#-------------------------------------------------------------------------------
# Replication code for survey 1

library(tidyverse); library(rio); library(cowplot); library(texreg)

# Set directory or change path to saved data
# Lines that save figures are commented out

# Data
# Respondents who did not consent already removed
tu <- import("survey1.csv")

#-------------------------------------------------------------------------------
# Initial cleaning

# Variable names to lower
colnames(tu) <- tolower(colnames(tu))

# Create and recode variables
tu <- tu %>%
  mutate(minutes = as.numeric(`duration (in seconds)`)/60, 
         pass_manip = ifelse(str_detect(abmanipluation, "Rising grocery"), 1, 0),
         believe = ifelse(str_detect(abbelievable, "unbeliev"), 0, 1)) %>%
  mutate(across(abeligible:abdoc2.1, ~ifelse(. == "", NA, .))) %>%
  
         # Treatment manipulation
  mutate(manipulation = case_when(rowSums(!is.na(across(abeligible:abmoney))) > 0  ~ "Baseline",
                                  rowSums(!is.na(across(abeligible.1:abdoc))) > 0  ~ "Min. Documentation",
                                  rowSums(!is.na(across(abeligible.2:abdoc.1))) > 0  ~ "Min. Documentation (est. Time)",
                                  rowSums(!is.na(across(abeligible.3:abdoc2))) > 0  ~ "Extensive Documentation",
                                  rowSums(!is.na(across(abeligible.4:abdoc2.1))) > 0  ~ "Extensive Documentation (est. Time)",
                                  T ~ NA_character_),
         # DVs
         # Believe eligible
         eligible_char = coalesce(abeligible, abeligible.1, abeligible.2, abeligible.3,
                                  abeligible.4),
         eligible = ifelse(eligible_char == "Yes", 1, 0),
         # Likely apply
         apply_char = coalesce(abapply_1, abapply_1.1, abapply_1.2, abapply_1.3, abapply_1.4),
         apply = case_when(apply_char == "Strongly disagree" ~ 1, 
                           apply_char == "Disagree" ~ 2, 
                           apply_char == "Somewhat disagree" ~ 3, 
                           apply_char == "Neither agree nor disagree" ~ 4,
                           apply_char == "Somewhat agree" ~ 5, 
                           apply_char == "Agree" ~ 6,
                           apply_char == "Strongly agree" ~ 7,
                           T ~ NA_real_),
         # Helpful for people like me
         helpful_char = coalesce(abapply_2, abapply_2.1, abapply_2.2, abapply_2.3, abapply_2.4),
         helpful = case_when(helpful_char == "Strongly disagree" ~ 1, 
                             helpful_char == "Disagree" ~ 2, 
                             helpful_char == "Somewhat disagree" ~ 3, 
                             helpful_char == "Neither agree nor disagree" ~ 4,
                             helpful_char == "Somewhat agree" ~ 5, 
                             helpful_char == "Agree" ~ 6,
                             helpful_char == "Strongly agree" ~ 7,
                             T ~ NA_real_),
         # Collecting information worth my time
         worth_char = coalesce(abapply_3, abapply_3.1, abapply_3.2, abapply_3.3, abapply_3.4),
         worth = case_when(worth_char == "Strongly disagree" ~ 1, 
                           worth_char == "Disagree" ~ 2, 
                           worth_char == "Somewhat disagree" ~ 3, 
                           worth_char == "Neither agree nor disagree" ~ 4,
                           worth_char == "Somewhat agree" ~ 5, 
                           worth_char == "Agree" ~ 6,
                           worth_char == "Strongly agree" ~ 7,
                           T ~ NA_real_),
         # Learning more about program
         learn_char = coalesce(abapply_4, abapply_4.1, abapply_4.2, abapply_4.3, abapply_4.4),
         learn = case_when(learn_char == "Strongly disagree" ~ 1, 
                           learn_char == "Disagree" ~ 2, 
                           learn_char == "Somewhat disagree" ~ 3, 
                           learn_char == "Neither agree nor disagree" ~ 4,
                           learn_char == "Somewhat agree" ~ 5, 
                           learn_char == "Agree" ~ 6,
                           learn_char == "Strongly agree" ~ 7,
                           T ~ NA_real_),
         # Dichotomize
         apply_di = ifelse(apply >= 5, 1, 0),
         helpful_di = ifelse(helpful >= 5, 1, 0),
         worth_di = ifelse(worth >= 5, 1, 0),
         learn_di = ifelse(learn >= 5, 1, 0),
         # Enough support
         enough_char = coalesce(abmoney, abmoney.1, abmoney.2, abmoney.3, abmoney.4), 
         not_enough = ifelse(str_detect(enough_char, "too little"), 1, 0),
         # Documentation
         documentation_min = coalesce(abdoc, abdoc.1),
         documentation_ext = coalesce(abdoc2, abdoc2.1),
         # Time estimation
         est_time = coalesce(abesttime_1_text, abesttime_1_text.1),
         no_est_time = coalesce(abesttime, abesttime.1),
         # Race
         race = case_when(ethnicity == 1 ~ "White",
                          ethnicity == 2 ~ "Black",
                          ethnicity == 3 ~ "American Indian",
                          ethnicity %in% seq(4, 14, 1) ~ "Asian/Pacific Islander",
                          ethnicity == 15 ~ "Other",
                          T ~ "Other"),
         race = ifelse(hispanic == 1, race, "Hispanic"),
         female = ifelse(gender == 2, "Female", "Male"),
         education = ifelse(education < 1, NA, education),
         hhi = ifelse(hhi == 25, NA, hhi)
         )

# Add additional indicators 
tu <- tu %>%
  mutate(manipulation_time = case_when(manipulation == "Baseline" ~ "Baseline",
                                         !is.na(abesttime_1_text.1) &
                                         manipulation == "Extensive Documentation (est. Time)" ~ "Extensive Doc. (gave hours)",
                                       is.na(abesttime_1_text.1) &
                                         manipulation == "Extensive Documentation (est. Time)" ~ "Extensive Doc. (not sure)",
                                       !is.na(abesttime_1_text) &
                                         manipulation == "Min. Documentation (est. Time)" ~ "Minimal Doc. (gave hours)",
                                       is.na(abesttime_1_text) &
                                         manipulation == "Min. Documentation (est. Time)" ~ "Minimal Doc. (not sure)"))
table(tu$manipulation_time)
table(tu$manipulation_time, tu$no_est_time)
table(tu$manipulation, tu$manipulation_time)

#-------------------------------------------------------------------------------
# Summary information 

# Timing
summary(tu$minutes)
median(tu$minutes)
table(tu$minutes>30)[2]
table(tu$minutes>15)[2]
table(tu$minutes<3)[2]

# Indicator for fewer than 3 minutes
tu <- tu %>%
  mutate(exclude_time = ifelse(minutes < 3, 1, 0))

# N per treatment
table(tu$manipulation)

# Manipulation check
table(tu$manipulation, tu$pass_manip)
sum(table(tu$manipulation, tu$pass_manip)[,2])/sum(table(tu$manipulation, tu$pass_manip))

# Believability
table(tu$believe)
prop.table(table(tu$manipulation, tu$believe), 1)

#------------------
# Means for each DV
tu %>%
  select(eligible, apply, helpful, worth, learn, not_enough) %>%
  summary()
tu %>%
  select(eligible, apply, helpful, worth, learn, not_enough) %>%
  summarize(across(everything(), ~ sd(.x, na.rm = TRUE)))

#-------------------------------------------------------------------------------
# Distributions

# Reshape to long format
long_tu <- tu %>%
  select(apply, helpful, worth, learn) %>%
  pivot_longer(cols = everything(), names_to = "question", values_to = "response") %>%
  filter(!is.na(response))

# Add labels
long_tu <- long_tu %>%
  mutate(
    question = recode(question,
                      apply = "Apply",
                      helpful = "Helpful for Me",
                      worth = "Worth Time",
                      learn = "Learn More"),
    response = factor(response, levels = 1:7,
                      labels = c("Strongly Disagree", "Disagree", "Somewhat Disagree",
                                 "Neither", "Somewhat Agree", "Agree", "Strongly Agree"))
  )

# Create the plot
ggplot(long_tu, aes(x = response)) +
  geom_bar(position = "dodge") +
  facet_wrap(~ question, ncol = 1) +
  labs(x = "", y = "Number of Respondents") +
  theme_bw(base_size = 12) +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1))
# ggsave("distribution_s1.pdf", width = 10, height = 8)
  
#-------------------------------------------------------------------------------
# Difference of means

# apply, helpful, worth, learn

# Function to grab the difference of means between two variables
dif.of.means.fun = function(var1, var2){
  temp <- t.test(var1, var2)
  p <- temp$p.value
  mean1 <- temp$estimate[1]
  mean2 <- temp$estimate[2]
  lower <- temp$conf.int[1]
  upper <- temp$conf.int[2]
  temp2 <- t.test(var1, var2,conf.level = 0.9)
  lower90 <- temp2$conf.int[1]
  upper90 <- temp2$conf.int[2]
  difference <- mean1 - mean2 
  out <- data.frame("mean1" = mean1,
                   "mean2"= mean2,
                   "difference" = difference,
                   "lower95" = lower,
                   "upper95" = upper,
                   "lower90" = lower90,
                   "upper90" = upper90,
                   "p" = p)
  return(out)
}

# Sanity check
dif.of.means.fun(tu$apply[tu$manipulation == "Min. Documentation"],
                 tu$apply[tu$manipulation == "Baseline"])
dif.of.means.fun(tu$apply[tu$manipulation == "Extensive Documentation (est. Time)"],
                 tu$apply[tu$manipulation == "Extensive Documentation"])

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
                             dif.of.means.fun(tu$apply[tu$manipulation == "Min. Documentation"],
                                              tu$apply[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$apply[tu$manipulation == "Min. Documentation (est. Time)"],
                                              tu$apply[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$apply[tu$manipulation == "Extensive Documentation"],
                                              tu$apply[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$apply[tu$manipulation == "Extensive Documentation (est. Time)"],
                                              tu$apply[tu$manipulation == "Baseline"]),
                             # helpful
                             dif.of.means.fun(tu$helpful[tu$manipulation == "Min. Documentation"],
                                              tu$helpful[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$helpful[tu$manipulation == "Min. Documentation (est. Time)"],
                                              tu$helpful[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$helpful[tu$manipulation == "Extensive Documentation"],
                                              tu$helpful[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$helpful[tu$manipulation == "Extensive Documentation (est. Time)"],
                                              tu$helpful[tu$manipulation == "Baseline"]),
                             # worth
                             dif.of.means.fun(tu$worth[tu$manipulation == "Min. Documentation"],
                                              tu$worth[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$worth[tu$manipulation == "Min. Documentation (est. Time)"],
                                              tu$worth[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$worth[tu$manipulation == "Extensive Documentation"],
                                              tu$worth[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$worth[tu$manipulation == "Extensive Documentation (est. Time)"],
                                              tu$worth[tu$manipulation == "Baseline"]),
                             # learn
                             dif.of.means.fun(tu$learn[tu$manipulation == "Min. Documentation"],
                                              tu$learn[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$learn[tu$manipulation == "Min. Documentation (est. Time)"],
                                              tu$learn[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$learn[tu$manipulation == "Extensive Documentation"],
                                              tu$learn[tu$manipulation == "Baseline"]),
                             dif.of.means.fun(tu$learn[tu$manipulation == "Extensive Documentation (est. Time)"],
                                              tu$learn[tu$manipulation == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. - \nBaseline",
                                             "Minimal Doc. (est. Time) - \nBaseline",
                                             "Extensive Doc. - \nBaseline",
                                             "Extensive Doc. (est. Time) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                          aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

#-------------------------------------------------------------------------------
# Difference of proportions

# eligible, not_enough

# Function to grab difference of proportions
prop.test.values <- function(generated.table, rounding.number = 3){
  n <- sum(generated.table)
  temp <- prop.test(generated.table)
  prop.1 <- temp$estimate[1]
  prop.2 <- temp$estimate[2]
  dif <- prop.1 - prop.2
  lower <- temp$conf.int[1]
  upper <- temp$conf.int[2]
  p <- temp$p.value
  temp2 <- prop.test(generated.table,conf.level = 0.9)
  lower90 <- temp2$conf.int[1]
  upper90 <- temp2$conf.int[2]
  out <- data.frame("n" = n,
                   "prop1" = prop.1,
                   "prop2" = prop.2,
                   "difference" = dif,
                   "lower95" = lower,
                   "upper95" = upper,
                   "lower90" = lower90,
                   "upper90" = upper90,
                   "p" = p)
  return(out)
}

# Sanity check
prop.test.values(table(tu$manipulation, tu$eligible)[c(1,2), ])

# Save as table
tab.elig <- table(tu$manipulation, tu$eligible)
tab.support <- table(tu$manipulation, tu$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Extensive Doc. - \nBaseline",
                       "Extensive Doc. (est. Time) - \nBaseline",
                       "Minimal Doc. - \nBaseline",
                       "Minimal Doc. (est. Time) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Minimal Doc. - \nBaseline",
                                 "Minimal Doc. (est. Time) - \nBaseline",
                                 "Extensive Doc. - \nBaseline",
                                 "Extensive Doc. (est. Time) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
       aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("figure1.png", width = 16, height = 8, dpi = 300)

#-------------------------------------------------------------------------------
# Difference of means
# Heterogeneous effects
# Split the estimation groups based on providing an estimate

# apply, helpful, worth, learn

# Sanity check
dif.of.means.fun(tu$apply[tu$manipulation_time == "Minimal Doc. (not sure)"],
                 tu$apply[tu$manipulation_time == "Baseline"])
dif.of.means.fun(tu$apply[tu$manipulation_time == "Minimal Doc. (gave hours)"],
                 tu$apply[tu$manipulation_time == "Baseline"])

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu$apply[tu$manipulation_time == "Minimal Doc. (not sure)"],
                   tu$apply[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$apply[tu$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu$apply[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$apply[tu$manipulation_time == "Extensive Doc. (not sure)"],
                   tu$apply[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$apply[tu$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu$apply[tu$manipulation_time == "Baseline"]),
  # helpful
  dif.of.means.fun(tu$helpful[tu$manipulation_time == "Minimal Doc. (not sure)"],
                   tu$helpful[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$helpful[tu$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu$helpful[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$helpful[tu$manipulation_time == "Extensive Doc. (not sure)"],
                   tu$helpful[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$helpful[tu$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu$helpful[tu$manipulation_time == "Baseline"]),
  # worth
  dif.of.means.fun(tu$worth[tu$manipulation_time == "Minimal Doc. (not sure)"],
                   tu$worth[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$worth[tu$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu$worth[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$worth[tu$manipulation_time == "Extensive Doc. (not sure)"],
                   tu$worth[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$worth[tu$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu$worth[tu$manipulation_time == "Baseline"]),
  # learn
  dif.of.means.fun(tu$learn[tu$manipulation_time == "Minimal Doc. (not sure)"],
                   tu$learn[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$learn[tu$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu$learn[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$learn[tu$manipulation_time == "Extensive Doc. (not sure)"],
                   tu$learn[tu$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu$learn[tu$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu$learn[tu$manipulation_time == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. (not sure) - \nBaseline",
                                             "Minimal Doc. (gave hours) - \nBaseline",
                                             "Extensive Doc. (not sure) - \nBaseline",
                                             "Extensive Doc. (gave hours) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

#-------------------------------------------------------------------------------
# Difference of proportions
# Heterogeneous effects

# eligible, not_enough

# Sanity check
prop.test.values(table(tu$manipulation_time, tu$eligible)[c(1,2), ])

# Save as table
tab.elig <- table(tu$manipulation_time, tu$eligible)
tab.support <- table(tu$manipulation_time, tu$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Ext. Doc. (gave hours) - \nBaseline",
                             "Ext. Doc. (not sure) - \nBaseline",
                             "Min. Doc. (gave hours) - \nBaseline",
                             "Min. Doc. (not sure) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Min. Doc. (not sure) - \nBaseline",
                                 "Min. Doc. (gave hours) - \nBaseline",
                                 "Ext. Doc. (not sure) - \nBaseline",
                                 "Ext. Doc. (gave hours) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("figure2.png", width = 17, height = 8, dpi = 300)

#-------------------------------------------------------------------------------
# Time estimation

# Just in estimating time category
tu.time <- tu %>%
  filter(str_detect(manipulation, "Time")) %>%
  mutate(manipulation2 = ifelse(str_detect(manipulation, "Min"), "Minimal Documention", "Extensive Documentation"),
         no_est_time = ifelse(no_est_time == "I'm not sure", 1, 0))

# Calculate the median and mean for each group
stats <- tu.time %>%
  filter(est_time < 11) %>%
  group_by(manipulation2) %>%
  summarize(
    median_est_time = median(est_time, na.rm = TRUE),
    mean_est_time = mean(est_time, na.rm = TRUE),
    sd_est_time = sd(est_time, na.rm   = TRUE)
  )

# Plot the histogram with a vertical line for the median and text labels for both median and mean
ggplot(filter(tu.time, est_time < 11), aes(est_time)) +
  geom_histogram(binwidth = 1, alpha = 0.5) +
  facet_wrap(~fct_rev(manipulation2)) +
  geom_vline(data = stats, aes(xintercept = median_est_time), color = "black", linetype = "dashed", linewidth = 1) +
  geom_text(
    data = stats, 
    aes(x = median_est_time, y = Inf, label = paste("Median:", round(median_est_time, 1))), 
    vjust = 5, hjust = -.5, color = "black", size = 4
  ) +
  geom_text(
    data = stats, 
    aes(x = median_est_time, y = Inf, label = paste("Mean:", round(mean_est_time, 1))), 
    vjust = 7, hjust = -.5, color = "darkblue", size = 4
  ) +
  theme_bw(base_size = 14) +
  labs(y = "", x = "Estimated Time")
# ggsave("figure3.png", width = 9, height = 4, dpi = 300)

#-------------------------------------------------------------------------------
# Documentation 

# Documentation variables
tu.min.doc <- tu %>%
  filter(manipulation %in% c("Min. Documentation", "Min. Documentation (est. Time)")) %>%
  mutate(equal = ifelse(documentation_min == "About equal", 1, 0),
         notsure = ifelse(documentation_min == "I'm not sure", 1, 0),
         identity = ifelse(documentation_min == "Proof of identity", 1, 0),
         income = ifelse(documentation_min == "Proof of income", 1, 0),
         benefits = ifelse(documentation_min == "Proof of other benefits", 1, 0),
         residency = ifelse(documentation_min == "Proof of residency", 1, 0))
tu.ext.doc <- tu %>%
  filter(manipulation %in% c("Extensive Documentation", "Extensive Documentation (est. Time)")) %>%
  mutate(equal = ifelse(documentation_ext == "About equal", 1, 0),
         notsure = ifelse(documentation_ext == "I'm not sure", 1, 0),
         identity = ifelse(documentation_ext == "Proof of identity", 1, 0),
         income = ifelse(documentation_ext == "Proof of income", 1, 0),
         benefits = ifelse(documentation_ext == "Proof of other benefits", 1, 0),
         residency = ifelse(documentation_ext == "Proof of residency", 1, 0),
         household = ifelse(documentation_ext == "Proof of household resources", 1, 0))


# Difference of proportion
# For each treatment group and outcome
diff.prop.min <- bind_rows(prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$equal)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$notsure)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$identity)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$income)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$benefits)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$residency)))
diff.prop.ext <- bind_rows(prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$equal)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$notsure)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$identity)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$income)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$benefits)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$residency)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$household)))

# Add in comparison text
diff.prop.min$compare <- "Estimated Time - \nMinimal Documentation"
diff.prop.ext$compare <- "Estimated Time - \nExtensive Documentation"
diff.prop.min$condition <- c("About equal", "I'm not sure", "Proof of identity", 
                             "Proof of income", "Proof of other benefits", 
                             "Proof of residency")
diff.prop.ext$condition <- c("About equal", "I'm not sure", "Proof of identity", 
                             "Proof of income", "Proof of other benefits", 
                             "Proof of residency", "Proof of household resources")

# Combine and graph
diff.prop.documentation <- rbind(diff.prop.min,
                                 diff.prop.ext)
# Change factor levels
diff.prop.documentation$compare <- fct_relevel(diff.prop.documentation$compare,
                                               "Estimated Time - \nMinimal Documentation",
                                               "Estimated Time - \nExtensive Documentation")
diff.prop.documentation$condition <- fct_relevel(diff.prop.documentation$condition,
                                                "I'm not sure", "About equal",
                                                "Proof of residency", "Proof of identity",
                                                "Proof of income", "Proof of other benefits", 
                                                "Proof of household resources")

doc.diff.graph <- ggplot(diff.prop.documentation, aes(x = fct_rev(condition), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~compare) +
  theme_bw(base_size = 14) +
  coord_flip() +
  labs(y = "Difference in Proportion", x = ""); doc.diff.graph

# Combine minimal and extensive documentation data
combined_doc <- bind_rows(
  tu.min.doc %>% 
    select(equal, notsure, identity, income, benefits, residency) %>% 
    summarise(across(everything(), ~ mean(.x, na.rm = TRUE))) %>%
    pivot_longer(cols = everything(), names_to = "doc", values_to = "percentage") %>%
    mutate(condition = "Minimal Documentation"),
  
  tu.ext.doc %>% 
    select(equal, notsure, identity, income, benefits, residency, household) %>%
    summarise(across(everything(), ~ mean(.x, na.rm = TRUE))) %>%
    pivot_longer(cols = everything(), names_to = "doc", values_to = "percentage") %>%
    mutate(condition = "Extensive Documenation")
)

# Convert to percentage
combined_doc$percentage <- combined_doc$percentage * 100

# Documentation type
combined_doc$doc <- c("About equal", "I'm not sure", "Proof of identity", 
                      "Proof of income", "Proof of other benefits", 
                      "Proof of residency", "About equal", "I'm not sure", "Proof of identity", 
                      "Proof of income", "Proof of other benefits", 
                      "Proof of residency", "Proof of household resources")
combined_doc$doc <- fct_relevel(combined_doc$doc,
                                "I'm not sure", "About equal",
                                "Proof of residency", "Proof of identity",
                                "Proof of income", "Proof of other benefits", 
                                "Proof of household resources")

doc.per.plot <- ggplot(combined_doc, aes(x = fct_rev(doc), y = percentage)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(y = "Percentage (%)", x = "") +
  theme_bw(base_size = 14) +
  scale_y_continuous(limits = c(0, 40), breaks = seq(0, 40, 10)) +
  facet_wrap(~fct_rev(condition)) +
  coord_flip(); doc.per.plot

# Combine
plot_grid(doc.per.plot, doc.diff.graph, nrow = 1,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1))

# ggsave("figure4.png", width = 14, height = 6, dpi = 300)

#-------------------------------------------------------------------------------
# Robustness checks
#-------------------------------------------------------------------------------

#---------------------------------
# Dichotomous DVs 
# apply, helpful, worth, learn
# Both full models and heterogeneous effects

# Full models

# Sanity check
prop.test.values(table(tu$manipulation, tu$apply_di)[c(1,2), ])

# Save as table
tab.apply <- table(tu$manipulation, tu$apply_di)
tab.helpful <- table(tu$manipulation, tu$helpful_di)
tab.worth <- table(tu$manipulation, tu$worth_di)
tab.learn <- table(tu$manipulation, tu$learn_di)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.apply[c(1, 2),]),
                       prop.test.values(tab.apply[c(1, 3),]),
                       prop.test.values(tab.apply[c(1, 4),]),
                       prop.test.values(tab.apply[c(1, 5),]),
                       prop.test.values(tab.helpful[c(1, 2),]),
                       prop.test.values(tab.helpful[c(1, 3),]),
                       prop.test.values(tab.helpful[c(1, 4),]),
                       prop.test.values(tab.helpful[c(1, 5),]),
                       prop.test.values(tab.worth[c(1, 2),]),
                       prop.test.values(tab.worth[c(1, 3),]),
                       prop.test.values(tab.worth[c(1, 4),]),
                       prop.test.values(tab.worth[c(1, 5),]),
                       prop.test.values(tab.learn[c(1, 2),]),
                       prop.test.values(tab.learn[c(1, 3),]),
                       prop.test.values(tab.learn[c(1, 4),]),
                       prop.test.values(tab.learn[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Extensive Doc. - \nBaseline",
                             "Extensive Doc. (est. Time) - \nBaseline",
                             "Minimal Doc. - \nBaseline",
                             "Minimal Doc. (est. Time) - \nBaseline"), 4))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Minimal Doc. - \nBaseline",
                                 "Minimal Doc. (est. Time) - \nBaseline",
                                 "Extensive Doc. - \nBaseline",
                                 "Extensive Doc. (est. Time) - \nBaseline")
diff.prop$condition <- c(rep("Apply", 4), rep("Helpful for Me", 4), 
                         rep("Worth Time", 4), rep("Learn More", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition, ncol = 2) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph
# ggsave("di_dv_s1.pdf", width = 14, height = 7)

# Heterogeneous effects

# Save as table
tab.apply <- table(tu$manipulation_time, tu$apply_di)
tab.helpful <- table(tu$manipulation_time, tu$helpful_di)
tab.worth <- table(tu$manipulation_time, tu$worth_di)
tab.learn <- table(tu$manipulation_time, tu$learn_di)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.apply[c(1, 2),]),
                       prop.test.values(tab.apply[c(1, 3),]),
                       prop.test.values(tab.apply[c(1, 4),]),
                       prop.test.values(tab.apply[c(1, 5),]),
                       prop.test.values(tab.helpful[c(1, 2),]),
                       prop.test.values(tab.helpful[c(1, 3),]),
                       prop.test.values(tab.helpful[c(1, 4),]),
                       prop.test.values(tab.helpful[c(1, 5),]),
                       prop.test.values(tab.worth[c(1, 2),]),
                       prop.test.values(tab.worth[c(1, 3),]),
                       prop.test.values(tab.worth[c(1, 4),]),
                       prop.test.values(tab.worth[c(1, 5),]),
                       prop.test.values(tab.learn[c(1, 2),]),
                       prop.test.values(tab.learn[c(1, 3),]),
                       prop.test.values(tab.learn[c(1, 4),]),
                       prop.test.values(tab.learn[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Ext. Doc. (gave hours) - \nBaseline",
                             "Ext. Doc. (not sure) - \nBaseline",
                             "Min. Doc. (gave hours) - \nBaseline",
                             "Min. Doc. (not sure) - \nBaseline"), 4))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Min. Doc. (not sure) - \nBaseline",
                                 "Min. Doc. (gave hours) - \nBaseline",
                                 "Ext. Doc. (not sure) - \nBaseline",
                                 "Ext. Doc. (gave hours) - \nBaseline")
diff.prop$condition <- c(rep("Apply", 4), rep("Helpful for Me", 4), 
                         rep("Worth Time", 4), rep("Learn More", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition, ncol = 2) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph
# ggsave("di_dv_het_s1.pdf", width = 14, height = 7)

#---------------------------------
# Dropping inattentive respondents
# 1) Failed manipulation check
# 2) Took fewer than 3 minutes to complete survey
tu.ia <- tu %>%
  filter(pass_manip == 1, 
         minutes > 3)

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation == "Min. Documentation"],
                   tu.ia$apply[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation == "Min. Documentation (est. Time)"],
                   tu.ia$apply[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation == "Extensive Documentation"],
                   tu.ia$apply[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation == "Extensive Documentation (est. Time)"],
                   tu.ia$apply[tu.ia$manipulation == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation == "Min. Documentation"],
                   tu.ia$helpful[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation == "Min. Documentation (est. Time)"],
                   tu.ia$helpful[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation == "Extensive Documentation"],
                   tu.ia$helpful[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation == "Extensive Documentation (est. Time)"],
                   tu.ia$helpful[tu.ia$manipulation == "Baseline"]),
  # worth
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation == "Min. Documentation"],
                   tu.ia$worth[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation == "Min. Documentation (est. Time)"],
                   tu.ia$worth[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation == "Extensive Documentation"],
                   tu.ia$worth[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation == "Extensive Documentation (est. Time)"],
                   tu.ia$worth[tu.ia$manipulation == "Baseline"]),
  # learn
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation == "Min. Documentation"],
                   tu.ia$learn[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation == "Min. Documentation (est. Time)"],
                   tu.ia$learn[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation == "Extensive Documentation"],
                   tu.ia$learn[tu.ia$manipulation == "Baseline"]),
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation == "Extensive Documentation (est. Time)"],
                   tu.ia$learn[tu.ia$manipulation == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. - \nBaseline",
                                             "Minimal Doc. (est. Time) - \nBaseline",
                                             "Extensive Doc. - \nBaseline",
                                             "Extensive Doc. (est. Time) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# Save as table
tab.elig <- table(tu.ia$manipulation, tu.ia$eligible)
tab.support <- table(tu.ia$manipulation, tu.ia$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Extensive Doc. - \nBaseline",
                             "Extensive Doc. (est. Time) - \nBaseline",
                             "Minimal Doc. - \nBaseline",
                             "Minimal Doc. (est. Time) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Minimal Doc. - \nBaseline",
                                 "Minimal Doc. (est. Time) - \nBaseline",
                                 "Extensive Doc. - \nBaseline",
                                 "Extensive Doc. (est. Time) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_attentive.pdf", width = 16, height = 8)

# Heterogeneous effects

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.ia$apply[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.ia$apply[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.ia$apply[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$apply[tu.ia$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.ia$apply[tu.ia$manipulation_time == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.ia$helpful[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.ia$helpful[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.ia$helpful[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$helpful[tu.ia$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.ia$helpful[tu.ia$manipulation_time == "Baseline"]),
  # worth
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.ia$worth[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.ia$worth[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.ia$worth[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$worth[tu.ia$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.ia$worth[tu.ia$manipulation_time == "Baseline"]),
  # learn
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.ia$learn[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.ia$learn[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.ia$learn[tu.ia$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.ia$learn[tu.ia$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.ia$learn[tu.ia$manipulation_time == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. (not sure) - \nBaseline",
                                             "Minimal Doc. (gave hours) - \nBaseline",
                                             "Extensive Doc. (not sure) - \nBaseline",
                                             "Extensive Doc. (gave hours) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# eligible, not_enough

# Sanity check
prop.test.values(table(tu.ia$manipulation_time, tu.ia$eligible)[c(1,2), ])

# Save as table
tab.elig <- table(tu.ia$manipulation_time, tu.ia$eligible)
tab.support <- table(tu.ia$manipulation_time, tu.ia$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Ext. Doc. (gave hours) - \nBaseline",
                             "Ext. Doc. (not sure) - \nBaseline",
                             "Min. Doc. (gave hours) - \nBaseline",
                             "Min. Doc. (not sure) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Min. Doc. (not sure) - \nBaseline",
                                 "Min. Doc. (gave hours) - \nBaseline",
                                 "Ext. Doc. (not sure) - \nBaseline",
                                 "Ext. Doc. (gave hours) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
  geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_het_time_attentive.pdf", width = 17, height = 8)

#-----------------------------------
# Prior contact with social welfare
# Documentation 
tu.personal <- tu %>%
  filter(str_detect(personal, "Medicaid|Unemployment|SNAP|Medicare|TANF"))
nrow(tu.personal)/nrow(tu)

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation == "Min. Documentation"],
                   tu.personal$apply[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation == "Min. Documentation (est. Time)"],
                   tu.personal$apply[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation == "Extensive Documentation"],
                   tu.personal$apply[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation == "Extensive Documentation (est. Time)"],
                   tu.personal$apply[tu.personal$manipulation == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation == "Min. Documentation"],
                   tu.personal$helpful[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation == "Min. Documentation (est. Time)"],
                   tu.personal$helpful[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation == "Extensive Documentation"],
                   tu.personal$helpful[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation == "Extensive Documentation (est. Time)"],
                   tu.personal$helpful[tu.personal$manipulation == "Baseline"]),
  # worth
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation == "Min. Documentation"],
                   tu.personal$worth[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation == "Min. Documentation (est. Time)"],
                   tu.personal$worth[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation == "Extensive Documentation"],
                   tu.personal$worth[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation == "Extensive Documentation (est. Time)"],
                   tu.personal$worth[tu.personal$manipulation == "Baseline"]),
  # learn
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation == "Min. Documentation"],
                   tu.personal$learn[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation == "Min. Documentation (est. Time)"],
                   tu.personal$learn[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation == "Extensive Documentation"],
                   tu.personal$learn[tu.personal$manipulation == "Baseline"]),
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation == "Extensive Documentation (est. Time)"],
                   tu.personal$learn[tu.personal$manipulation == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. - \nBaseline",
                                             "Minimal Doc. (est. Time) - \nBaseline",
                                             "Extensive Doc. - \nBaseline",
                                             "Extensive Doc. (est. Time) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# Save as table
tab.elig <- table(tu.personal$manipulation, tu.personal$eligible)
tab.support <- table(tu.personal$manipulation, tu.personal$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Extensive Doc. - \nBaseline",
                             "Extensive Doc. (est. Time) - \nBaseline",
                             "Minimal Doc. - \nBaseline",
                             "Minimal Doc. (est. Time) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Minimal Doc. - \nBaseline",
                                 "Minimal Doc. (est. Time) - \nBaseline",
                                 "Extensive Doc. - \nBaseline",
                                 "Extensive Doc. (est. Time) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_personal.pdf", width = 16, height = 8)

# Just in estimating time category
tu.time <- tu.personal %>%
  filter(str_detect(manipulation, "Time")) %>%
  mutate(manipulation2 = ifelse(str_detect(manipulation, "Min"), "Minimal Documention", "Extensive Documentation"))

table(tu.time$est_time)

# Calculate the median and mean for each group
stats <- tu.time %>%
  filter(est_time < 11) %>%
  group_by(manipulation2) %>%
  summarize(
    median_est_time = median(est_time, na.rm = TRUE),
    mean_est_time = mean(est_time, na.rm = TRUE),
    sd_est_time = sd(est_time, na.rm   = TRUE)
  )

# Plot the histogram with a vertical line for the median and text labels for both median and mean
ggplot(filter(tu.time, est_time < 11), aes(est_time)) +
  geom_histogram(binwidth = 1, alpha = 0.5) +
  facet_wrap(~fct_rev(manipulation2)) +
  geom_vline(data = stats, aes(xintercept = median_est_time), color = "black", linetype = "dashed", size = 1) +
  geom_text(
    data = stats, 
    aes(x = median_est_time, y = Inf, label = paste("Median:", round(median_est_time, 1))), 
    vjust = 5, hjust = -.5, color = "black", size = 4
  ) +
  geom_text(
    data = stats, 
    aes(x = median_est_time, y = Inf, label = paste("Mean:", round(mean_est_time, 1))), 
    vjust = 7, hjust = -.5, color = "darkblue", size = 4
  ) +
  theme_bw(base_size = 14) +
  labs(y = "", x = "Estimated Time")
# ggsave("time_hist_prior.pdf", width = 9, height = 4)

# Documentation 

# Documentation variables
tu.min.doc <- tu.personal %>%
  filter(manipulation %in% c("Min. Documentation", "Min. Documentation (est. Time)")) %>%
  mutate(equal = ifelse(documentation_min == "About equal", 1, 0),
         notsure = ifelse(documentation_min == "I'm not sure", 1, 0),
         identity = ifelse(documentation_min == "Proof of identity", 1, 0),
         income = ifelse(documentation_min == "Proof of income", 1, 0),
         benefits = ifelse(documentation_min == "Proof of other benefits", 1, 0),
         residency = ifelse(documentation_min == "Proof of residency", 1, 0))
tu.ext.doc <- tu.personal %>%
  filter(manipulation %in% c("Extensive Documentation", "Extensive Documentation (est. Time)")) %>%
  mutate(equal = ifelse(documentation_ext == "About equal", 1, 0),
         notsure = ifelse(documentation_ext == "I'm not sure", 1, 0),
         identity = ifelse(documentation_ext == "Proof of identity", 1, 0),
         income = ifelse(documentation_ext == "Proof of income", 1, 0),
         benefits = ifelse(documentation_ext == "Proof of other benefits", 1, 0),
         residency = ifelse(documentation_ext == "Proof of residency", 1, 0),
         household = ifelse(documentation_ext == "Proof of household resources", 1, 0))

# Difference of proportion
# For each treatment group and outcome
diff.prop.min <- bind_rows(prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$equal)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$notsure)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$identity)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$income)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$benefits)),
                           prop.test.values(table(tu.min.doc$manipulation, tu.min.doc$residency)))
diff.prop.ext <- bind_rows(prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$equal)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$notsure)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$identity)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$income)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$benefits)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$residency)),
                           prop.test.values(table(tu.ext.doc$manipulation, tu.ext.doc$household)))

# Add in comparison text
diff.prop.min$compare <- "Estimated Time - \nMinimal Documentation"
diff.prop.ext$compare <- "Estimated Time - \nExtensive Documentation"
diff.prop.min$condition <- c("About equal", "I'm not sure", "Proof of identity", 
                             "Proof of income", "Proof of other benefits", 
                             "Proof of residency")
diff.prop.ext$condition <- c("About equal", "I'm not sure", "Proof of identity", 
                             "Proof of income", "Proof of other benefits", 
                             "Proof of residency", "Proof of household resources")

# Combine and graph
diff.prop.documentation <- rbind(diff.prop.min,
                                 diff.prop.ext)
# Change factor levels
diff.prop.documentation$compare <- fct_relevel(diff.prop.documentation$compare,
                                               "Estimated Time - \nMinimal Documentation",
                                               "Estimated Time - \nExtensive Documentation")
diff.prop.documentation$condition <- fct_relevel(diff.prop.documentation$condition,
                                                 "I'm not sure", "About equal",
                                                 "Proof of residency", "Proof of identity",
                                                 "Proof of income", "Proof of other benefits", 
                                                 "Proof of household resources")

ggplot(diff.prop.documentation, aes(x = fct_rev(condition), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~compare) +
  theme_bw(base_size = 14) +
  coord_flip() +
  labs(y = "Difference in Proportion", x = "")
# ggsave("documentation_diff_prior.pdf", width = 10, height = 5)

# Heterogeneous effects

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.personal$apply[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.personal$apply[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.personal$apply[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$apply[tu.personal$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.personal$apply[tu.personal$manipulation_time == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.personal$helpful[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.personal$helpful[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.personal$helpful[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$helpful[tu.personal$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.personal$helpful[tu.personal$manipulation_time == "Baseline"]),
  # worth
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.personal$worth[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.personal$worth[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.personal$worth[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$worth[tu.personal$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.personal$worth[tu.personal$manipulation_time == "Baseline"]),
  # learn
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.personal$learn[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.personal$learn[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.personal$learn[tu.personal$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.personal$learn[tu.personal$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.personal$learn[tu.personal$manipulation_time == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. (not sure) - \nBaseline",
                                             "Minimal Doc. (gave hours) - \nBaseline",
                                             "Extensive Doc. (not sure) - \nBaseline",
                                             "Extensive Doc. (gave hours) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# eligible, not_enough

# Sanity check
prop.test.values(table(tu.personal$manipulation_time, tu.personal$eligible)[c(1,2), ])

# Save as table
tab.elig <- table(tu.personal$manipulation_time, tu.personal$eligible)
tab.support <- table(tu.personal$manipulation_time, tu.personal$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Ext. Doc. (gave hours) - \nBaseline",
                             "Ext. Doc. (not sure) - \nBaseline",
                             "Min. Doc. (gave hours) - \nBaseline",
                             "Min. Doc. (not sure) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Min. Doc. (not sure) - \nBaseline",
                                 "Min. Doc. (gave hours) - \nBaseline",
                                 "Ext. Doc. (not sure) - \nBaseline",
                                 "Ext. Doc. (gave hours) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_het_personal.pdf", width = 17, height = 8)


#--------------------------------
# Not eligible by reported income
# hhi: under 14 (14 = $75,000 to $79,999)
# hhi: under 20 (20 = $125,000 to $149,999)

# Assuming household income under $125,000
tu.125 <- tu %>%
  filter(hhi < 20)

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.125$apply[tu.125$manipulation == "Min. Documentation"],
                   tu.125$apply[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$apply[tu.125$manipulation == "Min. Documentation (est. Time)"],
                   tu.125$apply[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$apply[tu.125$manipulation == "Extensive Documentation"],
                   tu.125$apply[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$apply[tu.125$manipulation == "Extensive Documentation (est. Time)"],
                   tu.125$apply[tu.125$manipulation == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation == "Min. Documentation"],
                   tu.125$helpful[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation == "Min. Documentation (est. Time)"],
                   tu.125$helpful[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation == "Extensive Documentation"],
                   tu.125$helpful[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation == "Extensive Documentation (est. Time)"],
                   tu.125$helpful[tu.125$manipulation == "Baseline"]),
  # worth
  dif.of.means.fun(tu.125$worth[tu.125$manipulation == "Min. Documentation"],
                   tu.125$worth[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$worth[tu.125$manipulation == "Min. Documentation (est. Time)"],
                   tu.125$worth[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$worth[tu.125$manipulation == "Extensive Documentation"],
                   tu.125$worth[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$worth[tu.125$manipulation == "Extensive Documentation (est. Time)"],
                   tu.125$worth[tu.125$manipulation == "Baseline"]),
  # learn
  dif.of.means.fun(tu.125$learn[tu.125$manipulation == "Min. Documentation"],
                   tu.125$learn[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$learn[tu.125$manipulation == "Min. Documentation (est. Time)"],
                   tu.125$learn[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$learn[tu.125$manipulation == "Extensive Documentation"],
                   tu.125$learn[tu.125$manipulation == "Baseline"]),
  dif.of.means.fun(tu.125$learn[tu.125$manipulation == "Extensive Documentation (est. Time)"],
                   tu.125$learn[tu.125$manipulation == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. - \nBaseline",
                                             "Minimal Doc. (est. Time) - \nBaseline",
                                             "Extensive Doc. - \nBaseline",
                                             "Extensive Doc. (est. Time) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# Save as table
tab.elig <- table(tu.125$manipulation, tu.125$eligible)
tab.support <- table(tu.125$manipulation, tu.125$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Extensive Doc. - \nBaseline",
                             "Extensive Doc. (est. Time) - \nBaseline",
                             "Minimal Doc. - \nBaseline",
                             "Minimal Doc. (est. Time) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Minimal Doc. - \nBaseline",
                                 "Minimal Doc. (est. Time) - \nBaseline",
                                 "Extensive Doc. - \nBaseline",
                                 "Extensive Doc. (est. Time) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_125.pdf", width = 16, height = 8)

# Heterogeneous effects

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.125$apply[tu.125$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.125$apply[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$apply[tu.125$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.125$apply[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$apply[tu.125$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.125$apply[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$apply[tu.125$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.125$apply[tu.125$manipulation_time == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.125$helpful[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.125$helpful[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.125$helpful[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$helpful[tu.125$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.125$helpful[tu.125$manipulation_time == "Baseline"]),
  # worth
  dif.of.means.fun(tu.125$worth[tu.125$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.125$worth[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$worth[tu.125$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.125$worth[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$worth[tu.125$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.125$worth[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$worth[tu.125$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.125$worth[tu.125$manipulation_time == "Baseline"]),
  # learn
  dif.of.means.fun(tu.125$learn[tu.125$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.125$learn[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$learn[tu.125$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.125$learn[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$learn[tu.125$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.125$learn[tu.125$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.125$learn[tu.125$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.125$learn[tu.125$manipulation_time == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. (not sure) - \nBaseline",
                                             "Minimal Doc. (gave hours) - \nBaseline",
                                             "Extensive Doc. (not sure) - \nBaseline",
                                             "Extensive Doc. (gave hours) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# eligible, not_enough

# Sanity check
prop.test.values(table(tu.125$manipulation_time, tu.125$eligible)[c(1,2), ])

# Save as table
tab.elig <- table(tu.125$manipulation_time, tu.125$eligible)
tab.support <- table(tu.125$manipulation_time, tu.125$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Ext. Doc. (gave hours) - \nBaseline",
                             "Ext. Doc. (not sure) - \nBaseline",
                             "Min. Doc. (gave hours) - \nBaseline",
                             "Min. Doc. (not sure) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Min. Doc. (not sure) - \nBaseline",
                                 "Min. Doc. (gave hours) - \nBaseline",
                                 "Ext. Doc. (not sure) - \nBaseline",
                                 "Ext. Doc. (gave hours) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_het_time_125.pdf", width = 17, height = 8)

# Assuming income under $75,000
tu.75 <- tu %>%
  filter(hhi < 14)

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.75$apply[tu.75$manipulation == "Min. Documentation"],
                   tu.75$apply[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$apply[tu.75$manipulation == "Min. Documentation (est. Time)"],
                   tu.75$apply[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$apply[tu.75$manipulation == "Extensive Documentation"],
                   tu.75$apply[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$apply[tu.75$manipulation == "Extensive Documentation (est. Time)"],
                   tu.75$apply[tu.75$manipulation == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation == "Min. Documentation"],
                   tu.75$helpful[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation == "Min. Documentation (est. Time)"],
                   tu.75$helpful[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation == "Extensive Documentation"],
                   tu.75$helpful[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation == "Extensive Documentation (est. Time)"],
                   tu.75$helpful[tu.75$manipulation == "Baseline"]),
  # worth
  dif.of.means.fun(tu.75$worth[tu.75$manipulation == "Min. Documentation"],
                   tu.75$worth[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$worth[tu.75$manipulation == "Min. Documentation (est. Time)"],
                   tu.75$worth[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$worth[tu.75$manipulation == "Extensive Documentation"],
                   tu.75$worth[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$worth[tu.75$manipulation == "Extensive Documentation (est. Time)"],
                   tu.75$worth[tu.75$manipulation == "Baseline"]),
  # learn
  dif.of.means.fun(tu.75$learn[tu.75$manipulation == "Min. Documentation"],
                   tu.75$learn[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$learn[tu.75$manipulation == "Min. Documentation (est. Time)"],
                   tu.75$learn[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$learn[tu.75$manipulation == "Extensive Documentation"],
                   tu.75$learn[tu.75$manipulation == "Baseline"]),
  dif.of.means.fun(tu.75$learn[tu.75$manipulation == "Extensive Documentation (est. Time)"],
                   tu.75$learn[tu.75$manipulation == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. - \nBaseline",
                                             "Minimal Doc. (est. Time) - \nBaseline",
                                             "Extensive Doc. - \nBaseline",
                                             "Extensive Doc. (est. Time) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# Save as table
tab.elig <- table(tu.75$manipulation, tu.75$eligible)
tab.support <- table(tu.75$manipulation, tu.75$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Extensive Doc. - \nBaseline",
                             "Extensive Doc. (est. Time) - \nBaseline",
                             "Minimal Doc. - \nBaseline",
                             "Minimal Doc. (est. Time) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Minimal Doc. - \nBaseline",
                                 "Minimal Doc. (est. Time) - \nBaseline",
                                 "Extensive Doc. - \nBaseline",
                                 "Extensive Doc. (est. Time) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_75.pdf", width = 16, height = 8)

# Heterogeneous effects

# All treatments and DVs
dif.of.means.df <- bind_rows(# apply
  dif.of.means.fun(tu.75$apply[tu.75$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.75$apply[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$apply[tu.75$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.75$apply[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$apply[tu.75$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.75$apply[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$apply[tu.75$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.75$apply[tu.75$manipulation_time == "Baseline"]),
  # helpful
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.75$helpful[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.75$helpful[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.75$helpful[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$helpful[tu.75$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.75$helpful[tu.75$manipulation_time == "Baseline"]),
  # worth
  dif.of.means.fun(tu.75$worth[tu.75$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.75$worth[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$worth[tu.75$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.75$worth[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$worth[tu.75$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.75$worth[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$worth[tu.75$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.75$worth[tu.75$manipulation_time == "Baseline"]),
  # learn
  dif.of.means.fun(tu.75$learn[tu.75$manipulation_time == "Minimal Doc. (not sure)"],
                   tu.75$learn[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$learn[tu.75$manipulation_time == "Minimal Doc. (gave hours)"],
                   tu.75$learn[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$learn[tu.75$manipulation_time == "Extensive Doc. (not sure)"],
                   tu.75$learn[tu.75$manipulation_time == "Baseline"]),
  dif.of.means.fun(tu.75$learn[tu.75$manipulation_time == "Extensive Doc. (gave hours)"],
                   tu.75$learn[tu.75$manipulation_time == "Baseline"]))
# Outcome
dif.of.means.df$outcome <- c(rep(1, 4), rep(2, 4), rep(3, 4), rep(4, 4))
dif.of.means.df$outcome <- factor(dif.of.means.df$outcome,
                                  levels = 1:4, 
                                  labels = c("Apply", "Helpful for Me", "Worth Time", "Learn More"))

# Comparison
dif.of.means.df$compare <- c(rep(seq(1, 4), 4))
dif.of.means.df$compare <- factor(dif.of.means.df$compare,
                                  levels = 1:4, 
                                  labels = c("Minimal Doc. (not sure) - \nBaseline",
                                             "Minimal Doc. (gave hours) - \nBaseline",
                                             "Extensive Doc. (not sure) - \nBaseline",
                                             "Extensive Doc. (gave hours) - \nBaseline"))

# Plot
dif.of.mean.plot <- ggplot(dif.of.means.df, 
                           aes(x = fct_rev(compare), y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  theme_bw(base_size =  14)+
  facet_wrap(~outcome, ncol = 4) + 
  coord_flip() +
  labs(y = "Difference in Means", x = ""); dif.of.mean.plot

# eligible, not_enough

# Sanity check
prop.test.values(table(tu.75$manipulation_time, tu.75$eligible)[c(1,2), ])

# Save as table
tab.elig <- table(tu.75$manipulation_time, tu.75$eligible)
tab.support <- table(tu.75$manipulation_time, tu.75$not_enough)

# For each treatment group
diff.prop <- bind_rows(prop.test.values(tab.elig[c(1, 2),]),
                       prop.test.values(tab.elig[c(1, 3),]),
                       prop.test.values(tab.elig[c(1, 4),]),
                       prop.test.values(tab.elig[c(1, 5),]),
                       prop.test.values(tab.support[c(1, 2),]),
                       prop.test.values(tab.support[c(1, 3),]),
                       prop.test.values(tab.support[c(1, 4),]),
                       prop.test.values(tab.support[c(1, 5),]))

# Add in comparison text
diff.prop$compare <- c(rep(c("Ext. Doc. (gave hours) - \nBaseline",
                             "Ext. Doc. (not sure) - \nBaseline",
                             "Min. Doc. (gave hours) - \nBaseline",
                             "Min. Doc. (not sure) - \nBaseline"), 2))
diff.prop$compare <- fct_relevel(diff.prop$compare,
                                 "Min. Doc. (not sure) - \nBaseline",
                                 "Min. Doc. (gave hours) - \nBaseline",
                                 "Ext. Doc. (not sure) - \nBaseline",
                                 "Ext. Doc. (gave hours) - \nBaseline")
diff.prop$condition <- c(rep("Eligible for Program", 4), rep("Not Enough Support", 4))

# Graph
diff.prop.graph <- ggplot(diff.prop, 
                          aes(x = compare, y = difference)) +
  geom_hline(yintercept = 0, color = "black", linetype = "dashed")+
  geom_errorbar(width = 0, linewidth = 0.75, color = "darkgrey",
                aes(ymin = lower95, ymax = upper95)) +
 geom_errorbar(width = 0, linewidth=1.5,color="black",
                aes(ymin=lower90, ymax=upper90)) +
  geom_point(size = 5, fill = "black") +
  facet_wrap(~condition) +
  theme_bw(base_size = 14) +
  labs(y = "Difference in Proportion", x = ""); diff.prop.graph

# Combine graphs
plot_grid(diff.prop.graph, dif.of.mean.plot, nrow = 2,
          align = "hv", axis = "l", labels = c("a.", "b."),
          rel_heights = c(1, 1.2))
# ggsave("composite_main_het_time_75.pdf", width = 17, height = 8)

#-------------------------------------------------------------------------------
# Different moderators
# Perceive eligibility
# Adequate support
# Income

library(modelsummary)

# Low income indicator
# 5 is up to 34,999
tu <- tu %>%
  mutate(low_income = ifelse(hhi <= 5, 1, 0))
table(tu$low_income)

reg.elig.apply <- lm(apply ~ manipulation*eligible, data = tu); summary(reg.elig.apply)
reg.elig.worth <- lm(worth ~ manipulation*eligible, data = tu); summary(reg.elig.worth)
reg.elig.helpful <- lm(helpful ~ manipulation*eligible, data = tu); summary(reg.elig.helpful)
reg.elig.learn <- lm(learn ~ manipulation*eligible, data = tu); summary(reg.elig.learn)

texreg(
  l = list(reg.elig.apply, reg.elig.worth, reg.elig.helpful, reg.elig.learn),
  custom.model.names = c("Apply", "Helpful", "Worth Time", "Learn More"),
  stars = c(0.01, 0.05),
  digits = 2,
  booktabs = TRUE,
  use.packages = FALSE,
  caption = "Treatment moderated by perceived eligibility",
  label = "tab:mode_elig_s1"
)

reg.sup.apply <- lm(apply ~ manipulation*not_enough, data = tu); summary(reg.sup.apply)
reg.sup.worth <- lm(worth ~ manipulation*not_enough, data = tu); summary(reg.sup.worth)
reg.sup.helpful <- lm(helpful ~ manipulation*not_enough, data = tu); summary(reg.sup.helpful)
reg.sup.learn <- lm(learn ~ manipulation*not_enough, data = tu); summary(reg.sup.learn)

texreg(
  l = list(reg.sup.apply, reg.sup.worth, reg.sup.helpful, reg.sup.learn),
  custom.model.names = c("Apply", "Helpful", "Worth Time", "Learn More"),
  stars = c(0.01, 0.05),
  digits = 2,
  booktabs = TRUE,
  use.packages = FALSE,
  caption = "Treatment moderated by support perceived as not enough",
  label = "tab:mode_sup_s1"
)

reg.inc.apply <- lm(apply ~ manipulation*low_income, data = tu); summary(reg.inc.apply)
reg.inc.worth <- lm(worth ~ manipulation*low_income, data = tu); summary(reg.inc.worth)
reg.inc.helpful <- lm(helpful ~ manipulation*low_income, data = tu); summary(reg.inc.helpful)
reg.inc.learn <- lm(learn ~ manipulation*low_income, data = tu); summary(reg.inc.learn)

texreg(
  l = list(reg.inc.apply, reg.inc.worth, reg.inc.helpful, reg.inc.learn),
  custom.model.names = c("Apply", "Helpful", "Worth Time", "Learn More"),
  stars = c(0.01, 0.05),
  digits = 2,
  booktabs = TRUE,
  use.packages = FALSE,
  caption = "Treatment moderated by low-income ($< 35k$)",
  label = "tab:mode_inc_s1"
)

reg.elig.apply2 <- lm(apply ~ manipulation_time*eligible, data = tu); summary(reg.elig.apply)
reg.elig.worth2 <- lm(worth ~ manipulation_time*eligible, data = tu); summary(reg.elig.worth)
reg.elig.helpful2 <- lm(helpful ~ manipulation_time*eligible, data = tu); summary(reg.elig.helpful)
reg.elig.learn2 <- lm(learn ~ manipulation_time*eligible, data = tu); summary(reg.elig.learn)

texreg(
  l = list(reg.elig.apply2, reg.elig.worth2, reg.elig.helpful2, reg.elig.learn2),
  custom.model.names = c("Apply", "Helpful", "Worth Time", "Learn More"),
  stars = c(0.01, 0.05),
  digits = 2,
  booktabs = TRUE,
  use.packages = FALSE,
  caption = "Effects by Documentation Availability: Treatment moderated by perceived eligibility",
  label = "tab:mode_elig2_s1"
)

reg.sup.apply2 <- lm(apply ~ manipulation_time*not_enough, data = tu); summary(reg.sup.apply2)
reg.sup.worth2 <- lm(worth ~ manipulation_time*not_enough, data = tu); summary(reg.sup.worth2)
reg.sup.helpful2 <- lm(helpful ~ manipulation_time*not_enough, data = tu); summary(reg.sup.helpful2)
reg.sup.learn2 <- lm(learn ~ manipulation_time*not_enough, data = tu); summary(reg.sup.learn2)

texreg(
  l = list(reg.sup.apply2, reg.sup.worth2, reg.sup.helpful2, reg.sup.learn2),
  custom.model.names = c("Apply", "Helpful", "Worth Time", "Learn More"),
  stars = c(0.01, 0.05),
  digits = 2,
  booktabs = TRUE,
  use.packages = FALSE,
  caption = "Effects by Documentation Availability: Treatment moderated by support perceived as not enough",
  label = "tab:mode_sup2_s1"
)

reg.inc.apply2 <- lm(apply ~ manipulation_time*low_income, data = tu); summary(reg.inc.apply2)
reg.inc.worth2 <- lm(worth ~ manipulation_time*low_income, data = tu); summary(reg.inc.worth2)
reg.inc.helpful2 <- lm(helpful ~ manipulation_time*low_income, data = tu); summary(reg.inc.helpful2)
reg.inc.learn2 <- lm(learn ~ manipulation_time*low_income, data = tu); summary(reg.inc.learn2)

texreg(
  l = list(reg.inc.apply2, reg.inc.worth2, reg.inc.helpful2, reg.inc.learn2),
  custom.model.names = c("Apply", "Helpful", "Worth Time", "Learn More"),
  stars = c(0.01, 0.05),
  digits = 2,
  booktabs = TRUE,
  use.packages = FALSE,
  caption = "Effects by Documentation Availability: Treatment moderated by low-income ($< 35k$)",
  label = "tab:mode_inc2_s1"
)

#-------------------------------------------------------------------------------
# Balance checks

tu.balance <- tu %>%
  mutate(
    baseline = ifelse(manipulation == "Baseline", 1, 0),
    min_doc = ifelse(manipulation == "Min. Documentation", 1, 0),
    min_doc_est_time = ifelse(manipulation == "Min. Documentation (est. Time)", 1, 0),
    ext_doc = ifelse(manipulation == "Extensive Documentation", 1, 0),
    ext_doc_est_time = ifelse(manipulation == "Extensive Documentation (est. Time)", 1, 0)
  )


reg.balance1 = lm(baseline ~ factor(female) +
                    factor(race) +
                    as.numeric(age) + factor(party) + 
                    factor(urbanrural) + factor(region) +
                    as.numeric(hhi) + as.numeric(education), 
                  data = tu.balance)
reg.balance2 = lm(min_doc ~ factor(female) +
                    factor(race) +
                    as.numeric(age) + factor(party) + 
                    factor(urbanrural) + factor(region) +
                    as.numeric(hhi) + as.numeric(education), 
                  data = tu.balance)
reg.balance3 = lm(min_doc_est_time ~ factor(female) +
                    factor(race) +
                    as.numeric(age) + factor(party) + 
                    factor(urbanrural) + factor(region) +
                    as.numeric(hhi) + as.numeric(education), 
                  data = tu.balance)
reg.balance4 = lm(ext_doc ~ factor(female) +
                    factor(race) +
                    as.numeric(age) + factor(party) + 
                    factor(urbanrural) + factor(region) +
                    as.numeric(hhi) + as.numeric(education), 
                  data = tu.balance)
reg.balance5 = lm(ext_doc_est_time ~ factor(female) +
                    factor(race) +
                    as.numeric(age) + factor(party) + 
                    factor(urbanrural) + factor(region) +
                    as.numeric(hhi) + as.numeric(education), 
                  data = tu.balance)

texreg(list(reg.balance1, reg.balance2, reg.balance3, reg.balance4, reg.balance5))

#-------------------------------------------------------------------------------
# ANOVA and Regression

# DVs
dvs <- c("eligible", "apply", "helpful", "worth", "learn", "not_enough")

# Run ANOVAs for each DV and store results in a list
anova_results <- lapply(dvs, function(dv) {
  formula <- as.formula(paste(dv, "~ manipulation"))
  aov(formula, data = tu)
})

lapply(anova_results, summary)

# Run regressions for each dependent variable and store results in a list
regression_results <- lapply(dvs, function(dv) {
  formula <- as.formula(paste(dv, "~ min_doc + min_doc_est_time + 
                    ext_doc + ext_doc_est_time + 
                    factor(female) + factor(race) +
                    as.numeric(age) + factor(party) + 
                    factor(urbanrural) + factor(region) +
                    as.numeric(hhi) + as.numeric(education)"))
  lm(formula, data = tu.balance)
})

r.results <- lapply(regression_results, summary)
texreg(r.results)

#-------------------------------------------------------------------------------
# Summary statistics
tu.balance %>%
  select(baseline, min_doc, min_doc_est_time, ext_doc, ext_doc_est_time, 
         eligible, apply, helpful, worth, learn, not_enough,
         female, race, age, party, urbanrural, region, hhi, education, 
         pass_manip, minutes) %>%
  stargazer::stargazer(type = "latex", summary = T, digits = 2)
